function fv=plm(y,xs,r0,a,b)
%******************************************************************
%** This is the core program of truncated regression model with  **
%** constrained optimization (TRMCO)                             **
%** Created by XXXXXXXXXXX, Dec 19, 2011                         **
%** specifically written for Simulation II                       ** 
%** y:dependent variable (n by 1)                                **
%** xs:covariate matrix without the constant (n by m)            **
%** r0:initial value of parameters (m+2 by 1)                    **
%** a:lower boundary limit (scalar)                              **
%** b:upper boundary limit (scalar)                              **
%** output:parameter estimates, negative inverse Hessian,        **
%** and related statistics                                       **  
%******************************************************************
format long;
% basic information
warning off all; n=size(y,1); stop=1; 
m=size(xs,2); record=ones(1,4*m+17); 
% initial values
lamda=ones(2*m+6,1)*10; kk=1; run=1; runmax=101; x=[ones(n,1) xs];
%****************** Record the initial results *******************************
sumD1=0;
  for i=1:1:n
        if abs(b-x(i,:)*r0(1:m+1))<5*r0(m+2) || abs(a-x(i,:)*r0(1:m+1))<5*r0(m+2)
           D=normcdf((b-x(i,:)*r0(1:m+1))/r0(m+2))-normcdf((a-x(i,:)*r0(1:m+1))/r0(m+2));           
        else
           D=1; 
        end   
  sumD1=sumD1-log(D*sqrt(2*pi)*r0(m+2));
  end   
record(1,1)=0;
record(1,2:m+3)=r0';
record(1,m+4:2*m+5)=zeros(1,4);
record(1,2*m+6:4*m+11)=lamda';
record(1,4*m+12)=sumD1-1/(2*r0(m+2)^2)*((y-x*r0(1:m+1))'*(y-x*r0(1:m+1)));
record(1,4*m+13:4*m+17)=ones(1,5)*-999;
% xmax and xmin
xmax=zeros(m+1,1); xmin=zeros(m+1,1);
   for i=1:1:m+1
       xmax(i)=max(x(:,i));
       xmin(i)=min(x(:,i));
   end    
xlimit=[xmax xmin];   
while stop>0 && run<runmax
% vp and vn
vp=zeros(m+1,1); vn=zeros(m+1,1);
for i=1:1:m+1
    if r0(i)>=0
       vp(i)=1; 
    else
       vn(i)=1; 
    end
end    
v=[vp vn];
% compute necessary matrices
ci=cpar(r0,a,b,x,xlimit,v); df=dfpar(r0,a,b,y,x); 
ai=aipar(r0,a,b,x,xlimit,v); 
%****************** Step 1:check KKT conditions ********************
% KKT
KKT1=df+ai'*lamda; check1=sum(abs(KKT1)>10^-4);
KKT2=ci; check2=sum(KKT2>10^-4);
KKT3=lamda; check3=sum(KKT3<-10^-4);
KKT4=lamda'*ci; check4=sum(abs(KKT4)>10^-4);
check=check1+check2+check3+check4;
%check=check1+check2+check4;
   if check==0
      stop=0; 
   end
% if not stoped, then continue iterating
%****************** Step 2:compute dk and lamdaQP ********************
L=Lpar(r0,a,b,y,x);   
% solving the QP problem  
d=quadprog(L,df,ai,-ci);  
lamdaQP=(ai')\(-df-L*d);
%****************** Step 3:set new x and lamda ***********************
r1=r0+d;
lamda=lamdaQP;
%****************** Record the results *******************************
sumD2=0;
  for i=1:1:n
        if abs(b-x(i,:)*r1(1:m+1))<5*r1(m+2) || abs(a-x(i,:)*r1(1:m+1))<5*r1(m+2)
           D=normcdf((b-x(i,:)*r1(1:m+1))/r1(m+2))-normcdf((a-x(i,:)*r1(1:m+1))/r1(m+2));           
        else
           D=1; 
        end   
  sumD2=sumD2-log(D*sqrt(2*pi)*r1(m+2));
  end  
record(kk+1,1)=kk;
record(kk+1,2:m+3)=r1';
record(kk+1,m+4:2*m+5)=d';
record(kk+1,2*m+6:4*m+11)=lamdaQP';
record(kk+1,4*m+12)=sumD2-1/(2*r1(m+2)^2)*((y-x*r1(1:m+1))'*(y-x*r1(1:m+1)));
record(kk+1,4*m+13)=record(kk+1,4*m+12)-record(kk,4*m+12);
record(kk+1,4*m+14:4*m+17)=[check1 check2 check3 check4];
%****************** Step 4 and 5: return to Step 1 *******************
r0=r1;
kk=kk+1;
run=run+1;
end
% check convergence
if run==runmax  
   maxLL=-10^20; maxt=-999;
   for t=2:1:runmax-1
        if record(t,4*m+12)>maxLL && record(t+1,4*m+15)==0
           maxLL=record(t,4*m+12);
           maxt=t;
        end
   end
  fs=record(maxt,2:m+3);
  note=record(maxt+1,4*m+15);
else
  if record(kk-1,m+3)<10^10
     fs=record(kk-1,2:m+3);   
     note=record(kk,4*m+15);
     maxt=kk-1;
  else
     maxLL=-10^20; maxt=-999; 
       for t=2:1:run-2
           if record(t,4*m+12)>maxLL && record(t+1,4*m+15)==0
              maxLL=record(t,4*m+12);
              maxt=t;
           end
       end
       if maxt~=-999 
          fs=record(maxt,2:m+3);
          note=record(maxt+1,4*m+15);
       end   
  end       
end
if maxt==-999
    fv=ones(6,4)*-999; 
else
    var=Lpar(fs',a,b,y,x)\eye(4);
    last=[run record(2,4*m+15) note 0];
    fv=[fs;var;last];
end    



